Self-organized criticality in a model of biological evolution with long range interactions 
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In this work we study the effects of introducing long range interactions in the Bak-Sneppen (BS) 
model of biological evolution. We analyze a recently proposed version of the BS model where the 
interactions decay as r~ a ; in this way the first nearest neighbors model is recovered in the limit 
a — > oo, and the random neighbors version for a = 0. We study the space and time correlations 
and analize how the critical behavior of the system changes with the parameter a. We also study 
the sensibility to initial conditions of the model using the spreading of damage technique. We found 
that the system displays several distinct critical regimes as a is varied from a = to a — » oo 



In recent years an increasing numbers of systems that 
present Self Organized Criticality have been widely 
investigated. The general approach of statistical physics, 
where simple models try to catch the essencial ingre- 
dientes responsable for a given complex behavior has 
turned out to be very powerful for the study of this kind 
of problems. In particular Bak and Sneppen || have 
introduced a simple model which has shown to be able 
to reproduce evolutionary features such as punctuacted 
equilibrium Altough this model does not intend to 
give an accurate description of darwinian evolution, it 
catches into a single and very simple scheme (it is based 
on very simple dynamical rules) several features that are 
expected to be present in evolutionary processes, that 
is, punctuated equilibrium fl, Self Organized Critical- 
ity (SOC) J3| and weak sensitivity to initial conditions 
(WSIC) l^ll (i.e., chaotic behaviour where the trajecto- 
ries depart with a power law of the time instead of expo- 
nentially). In this sense, one important question arises 
about the robustness of such properties against modifica- 
tions (i.e., complexifying) of the simple dynamical rules 
on which the model is based. The original model, here- 
after referred as the first nearest neighbors (FNN) ver- 
sion ||, includes only nearest neighbors interactions in 
a one dimensional chain. This model presents SOC ||] 
and weak sensitivity to initial conditions (^,^). On the 
other hand, another version of the model with interac- 
tions between sites randomly chosen in the lattice (and 
therefore it can be regarded as a mean field version of the 
FNN), hereafter referred as the random negihbors (RN) 
version Q], does not present SOC ||. Moreover, it is not 
expected (and we shall show in this work that it is indeed 
the case) to present WSIC. 

Systems of coevolutionary species are expected to have 
some distance decaying interactions, thus lying somehow 
between the two previous schemes. Although not well de- 
fined, the concept of "distance" between species in these 
scenarios may be regarded as associated to some com- 
plex network of relationships including competition for 
resources and predator-prey ones, among many others 



[|9|. In this sense, the environmental modifications pro- 
duced by the extinction of one species may be expected to 
affect many others not directly related to it, where the in- 
tensity of such influence depends on the above mentioned 
distance. 

Along this line, in this letter we will focuse on the ro- 
bustness of the SOC and sensitivity to initial conditions 
properties of the Bak and Sneppen model against the in- 
troduction of long-range distance dependent interactions. 
To this end, we consider a generalization of the model, 
recently proposed by Cafiero et al |l0| that takes into ac- 
count long-range interactions between species that decay 
as r~ a , where r represent the distance between species 
(mesured in lattice units, i.e., r = 1, 2, ... in a chain) and 
a > is a parameter that control s the effective range of 
the interactions. The major value of this gen eralization, 
unlike others introduced in the literature [ [1 l|Jl 2|| , resides 
in the fact that it allows simply to retrieve the two above 
mentioned models by varying continuously the parame- 
ter a: when a — > we recover the RN version while for 
a — ► oo we recover the FNN one. 

The model consist of an iV-site linear lattice with peri- 
odic boundary conditions (i.e., a ring of N sites), where 
each site represents a species. Each species has associated 
a real variable bj, < bj < 1, that measures the relative 
fitness barrier. Starting from a random barrier distribu- 
tion, at each successive time step we identify the smallest 
barrier bj , and modify it by choosing a new random value 
from a uniform distribution. This change represents a 
jump of a species across its fitness barrier to a mutated 
species. This mutation must also affect other species in 
the chain, and to take into account this phenomena one 
defines a neighborhood which will also be modified in 
the same way. In Ref. || (a — > oo) the authors consid- 
ered the case in which this neighborhood consists of the 
two nearest species of the mutating one while in Ref. Q| 
(a = 0) the neighborhood consist of if — 1 species chosen 
at random among all the species of the chain. 

In order to generalize these models, we choose the 
neighborhood (of K — 1 species) at random with a prob- 
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ability that decays as r^ Q , where r\j is the distance of a 
given species j to the species with the smallest barrier i, 
and giving them new random values chosen from a uni- 
form distribution. In this wayfor K = 3 and a^oowe 
recover the two nearest neighbors model, while for a = 
we reproduce the random neighbors mean field version. 

To determine whether the system attains a self- 
organized critical state, we analyze the following quan- 
tities: the barrier distribution P(b) in the final steady 
state, the spatial correlation C(r) and the first return 
time distribution C(t). In order to study the sensitiv- 
ity to initial conditions, we calculate the time evolu- 
tion of the Hamming distance D(t) between two different 
replicas submitted to the same noise (damage spreading 
method). 

Since our main interest is to analyze the crossover be- 
tween the limits a — and a — > oo, we will restrict our- 
selves to consider the K — 3 case. Latter on we will dis- 
cuss briefly the effect of increasing K. Figure [l] presents 
the distribution P{b) of barrier values, for three different 
a values. Note that, independently of a the curves are 
qualitatively similar. The typical behavior of these curves 
can be characterized by the value of P(b(l)) (i.e, the sat- 
uration value of the distribution), which is displayed in 
Fig.^ as a function of a for three different system sizes. 
We can clearly distinguish three different regimes. For 
a < 1 the value of P(b(l)) is independent of N and the 
behavior observed for P(b) colapses to the one observed 
in the RN Bak-Sneppen model. There exists some value 
a — a c such that for intermediate values 1 < a < a c the 
value of P(6(l)) is very sensitive to changes in a, increas- 
ing its value as a grows, and finally, for a > a c , P(b(l)) 
reaches a saturation value, and we recover the behavior 
of P{b) for the FNN model when TV — > oo. The value of 
a c can be roughly estimated from numerical extrapola- 
tions of the curves to 1/N — > 0. We obtained that a c « 4 
for K — 3 (further analysis of the critical exponents will 
confirm this estimation). 



FIG. 1. Distribution of barrier values for three different 
values of a : 0.5, 3.0, and 100.0. All cases have N = 250. 




FIG. 2. P(b(l)) vs a for three different system sizes 
N = 250, 500, and 1000. 

We consider now the spatial and temporal correlations 
between the minimum barriers in order to determine the 
presence of SOC. Figure |^ presents in a log-log plot the 
probability C(r) that the minimum barriers at two suc- 
cesive updates will be separated by r sites. We observed 
a power law behavior C(r) ~ r _7r for all a =/= 0. In Fig. |] 
we present how n changes with a. When a = the spa- 
tial correlation is constant (w = 0) as in the RN model. 
As a grows n increases until it reaches a saturation value 
7r = 3.2 ± 0.2 for a > a c , in agreement with the results 
observed in the FNN model. 
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FIG. 3. Log-log plot of the space correlation C(r), for four 
different values of a : 0.0, 0.5, 2, 5, and 10.0. The system size 
is N = 1000 
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FIG. 4. Space correlation exponent ir vs a. 

Next we calculate the first return time distribution 
C(t), defined, as the probability that, if a given site is 
the minimum at time to, it will again be the minimum 
for the first time at time to + t. In Fig. [| we present 
our results for four different values of a (0.5, 1.5, 2.0 and 
3.0). For a > 2 the first return time clearly presents a 
power law behavior C(t) ~ i -r , even for finite system 
sizes. 

For a < 2 the system displays finite size effects, as can 
be seen in Fig. |] where we present C(i) when a = 0.5 for 
different system sizes; it is clear that a power law decay 
emerges as N grows. 



FIG. 6. Log-log plot of the first return probability C(t), for 
a = 0.5, for three different sizes, N = 1000, 1500, and 2000. 

In Fig. ^ we show how the first return time exponent r 
depends on a. Here again we find three different regimes: 
For a < 1 (unlike the spatial correlation exponent) all the 
curves C a {t) colapse and r = 1.5, displaying the same be- 
havior found in the RN Bak-Sneppen model with K = 3, 
where r = 3/2 exactly For 1 < a < a c , the value of 
r strongly depends on the value of a, with a minimum 
for a ~ 1.6, in agreement with the results of Cafiero et al 
|]l0f . For a > a c the value of r attains a saturation value 
t = 1.56 ± 0.05 in agreement with the value observed in 
the FNN Bak-Sneppen model where r = 1.6. 
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FIG. 7. First return probability exponent, r, as a function 
of a. 

FIG. 5. Log-log plot of the first return probability C(t), for Summarizing the results displayed in these figures, for 

four different values of a : 0.5, 1.5,2.0, and 3.0. System size < a < 1, since r presents the same trivial value ob- 
is iV = 1000 served in the RN Bak-Sneppen model, we cannot regard 
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the system as critical [||. For 1 < a < a c the expo- 
nents depend strongly on a, and since the exponents are 
non trivial, we regard this as a strong indicator of crit- 
icality in the system. Finally for a > a c the exponents 
becomes independent of a, taking the short range values 
observed in the FNN Bak-Sneppen model. We have ob- 
served that as we increase the number of interacting sites 
K , the value of a c decreases, slowly converging to a c = 2. 
This behavior reminds that of the one dimensional ferro- 
magnetic Ising model with the same type of interactions 
presented here, where the borderline between short and 
long range critical regimes is a = 2 jl3j . 

Next, we study the sensibility to initial conditions of 
this model and its dependence on a. To do that, we use 
the spreading of damage technique, which had previouly 
been applied to the FNN model |^J. In this particular 
limit it was shown that the system presents a weak sen- 
sivity to initial conditions, characterized by a power law 
increment, as times goes on, of the Hamming distance 
between replicas of the system. This behavior is reminis- 
cent of those observed at the edge of chaos in dynamical 
systems with few degrees of freedom. The procedure is as 
follows: given a configuration of N barrier values ({&•• }) 
in the self-organized critical state, we create a replica of 

(2) 

the system ({bj }) by choosing a site randomly and in- 
terchanging the value of this site with the value of the 
site with the smallest barrier. From then on we use the 
same random numbers for updating the barrier values in 
both replicas. 

Wc define the Hamming distance between the two 
replicas as: 

1 N 

D(t) = -J2\ b ? ] (t)-bf\t)\ (1) 

If the Hamming distance goes to zero we say that the 
system is in a frozen phase. On the other hand if the 
Hamming distance remains non zero we say that the sys- 
tem is chaotic in analogy with dynamical systems. 

Regarding the behavior of the average normalized 
Hamming distance D(N, t) = (D(t)} / we ob- 

served two different regimes as a is varied. For a < 2, 
D(N,t) reaches a saturation value D(N, oo) in just one 
step. The quotient D(N, oo)/N — > when N — > oo, and 
therefore the system does not present sensibility to the 
initial conditions. 

The typical temporal behaviour of D(N,t) for a > 2 
is displayed in Fig.|| for a = 2.5 and three different sys- 
tem sizes (the results presented correspond to averages 
over 500 * N realizations). We see that D(N,t) ~ t 5 
when t <C N and it saturates into a system size depen- 
dent value for large times, clearly showing WSIC pLp|. 
Moreover, we verified for Va > 2 the finite size scaling 
behavior 



D(N,t) ~ N« F (J^j ( 2 ) 

where 5 = 0.40 @, and z = 1.7 ± 0.2 is the dynami- 
cal exponent defined by t s ~ N z , t s being the value of 
t at which the increasing regime crosses over onto the 
saturation regime pl| (given by the intersection of the 
linearly increasing branch of the curve and the horizon- 
tal branch). Both exponents are independent of a. 
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FIG. 8. Normalized Hamming distance D(N,t) for a — 2.5 
and three system sizes, N = 1000, 1500, and 2000. 

Concluding, we have studied how long range interac- 
tions affects the criticality of the stationary state of our 
model and its sensibility to the initial conditions. Con- 
cerning the SOC, we observed three different regimes de- 
pending on a. For a > a c we can speak about a short- 
range critical regime, where the system presents SOC. 
Moreover, we observe that this property displays uni- 
versality, in the sense that most of the associated crit- 
ical exponents are independent of a. For < a < 1 
the system does not present SOC, although C(r) dis- 
plays a non-trivial power law decay with r, unlike the 
RN model for wich C(r) is constant. Moreover, all the 
relevant state functions or distributions become indepen- 
dent of a. This behavior has already been observed in 
a variety of systems with long-range interactions, both 
related to equilibrium |l^ , [l7| and non-equilibrium prop- 
erties Jl8| , p^[ . In all these systems it has been observed 
that the mean field behavior becomes dominant when 
< a < d, d being the dimensionality of the underlying 
lattice. Hence, in our case we can speak about a "mean- 
field" (non-critical) regime, i.e., that of the RN model. 
Finally, for 1 < a < a c we have a long-range critical 
regime, where the system presents non-universal SOC, 
i.e., the associated critical exponents depend strongly on 
a. 
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Concerning the sensibility to the initial conditions we 
observed two regimes: one for < a < 2 where the sys- 
tem does not present sensibility to the initial conditions 
of any type, while for a > 2 it displays universal WSIC, 
in the sense that the exponents of the scaling law (Q) are 
independent of a. 

We see that, at least one of the borderline values (and 
probably all of them) that separate the different regimes 
seems to be directly related to the dimensionality of the 
system. Hence, such dimensionality appears as a fun- 
damental parameter to determine the robustness of the 
model against variations in the range of the interactions. 
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